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Abstract 

Starting from a factorization theorem in Soft-Collinear Effective Theory, the thrust dis- 
tribution in e + e~ collisions is calculated including resummation of the next-to-next-to- 
next-to leading logarithms. This is a significant improvement over previous calculations 
which were only valid to next-to-leading logarithmic order. The fixed-order expansion of 
the resummed result approaches the exact fixed-order distribution towards the kinematic 
endpoint. This close agreement provides a verification of both the effective field theory 
expression and recently completed next-to-next-to-leading fixed-order event shapes. The 
resummed distribution is then matched to fixed order, resulting in a distribution valid 
over a large range of thrust. A fit to aleph and OPAL data from lep 1 and LEP 2 
produces a s {m z ) = 0.1172 ± 0.0010 ± 0.0008 ± 0.0012 ± 0.0012, where the uncertainties 
are respectively statistical, systematic, hadronic, and perturbative. This is one of the 
world's most precise extractions of a s to date. 



1 Introduction 



Lepton colliders, such as the Large Electron-Positron collider lep which ran from 1989-2000 
at CERN, provide an optimal environment for precision studies in high energy physics. Lacking 
the complications of strongly interacting initial states, which plague hadron colliders, lep has 
been able to provide extremely accurate measurements of standard model quantities such as 
the Z-boson mass, and its results tightly constrain beyond-the-standard model physics. The 
precision lep data is also used for QCD studies, for example to determine the strong coupling 
constant a s . With the variation of a s known to 4- loops, one should be able to confirm in 
great detail the running of the coupling, or use it to establish a discrepancy which might 
indicate new physics. Even at fixed center-of-mass energy, differential distributions for event 
shapes, such as thrust probe several energy scales and are extremely sensitive to the running 
coupling. Moreover, event shape variables are designed to be infrared safe, so that they can be 
calculated in perturbation theory and so the theoretical predictions should be correspondingly 
clean. Nevertheless, extractions of a s from event shapes at lep have until now been limited 
by theoretical uncertainty from unknown higher order terms in the perturbative expansion. 

One difficulty in achieving an accurate theoretical prediction from QCD has been the 
complexity of the relevant fixed-order calculations. Indeed, while the next-to-leading-order 
(NLO) results for event shapes have been known since 1980 [lj, the relevant next-to- next- 
to-leading order (NNLO) calculations were completed only in 2007 [HE]- In addition to the 
loop integrals, the subtraction of soft and collinear divergencies in the real emission diagrams 
presented a major complication. In fact, this is the first calculation where a subtraction scheme 
has been successfully implemented at NNLO [I] . However, even with these new results at hand, 
the corresponding extraction of a s continues to be limited by perturbative uncertainty. The 
result of [5] was a s (mz) = 0.1240 ± 0.0033, with a perturbative uncertainty of 0.0029. This 
NNLO result for the strong coupling constant comes out lower than at NLO, but 2a higher 
than the PDG average a s (mz) = 0.1176 ± 0.0020 [6]. Actually, the most precise values of a s 
are currently determined not from lep but at low energies using lattice simulations [7] and 
r-decays [8]. An extensive review of a s determinations is given in j9], new determinations 
since its publication include [TOj [11] . 

To further reduce the theoretical uncertainty of event shape calculations, it is important 
to resum the dominant perturbative contributions to all orders in a s . To see this, consider 
thrust, which is defined as 

T = max ' — — , (1) 

where the sum is over all momentum 3- vectors p; in the event, and the maximum is over all 
unit 3-vectors n. In the endpoint region, T — > 1 or r = (1 — T) — > 0, no fixed-order calculation 
could possibly describe the full distribution due to the appearance of large logarithms. For 
example, at leading order in perturbation theory the thrust distribution has the form 

1 da . 2a s 
a dr Sty 

where the ellipsis denotes terms that are regular in the limit r — > 0. Upon integration over 
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the endpoint region, one finds 



i?(r)= r dr '!^ = l + ^[-21n 2 T-31nr + ...] . (3) 
Jo ffodr' 3tt 

Double logarithmic terms of the form a™ ln 2n r arise from regions of phase space where the 
quarks or gluons are soft or collinear. For small enough r, higher order terms are just as 
important as lower order ones and the standard perturbative expansion breaks down. Re- 
summation refers to summing a series of contributions of the form a™ ln m r for the integral 
R(r) or a™(ln m_1 r)/r for the differential distribution. Leading logarithmic (LL) accuracy is 
achieved by summing the tower of logarithms with m = 2n, next-to-leading logarithmic accu- 
racy (NLL) also sums the terms with m — 2n — 1. Resummation at N fc LL accuracy, provides 
all logarithmic terms with 2n > m > 2n — 2k + 1, as detailed in Section [2J 

The first resummation of event shapes was done by Catani, Trentadue, Turnock and Web- 
ber (CTTW) in [12J. Their approach was to define jet functions Jc(p 2 ) as the probability 
for finding a jet of invariant mass p 2 in the event. These can be calculated to NLL by sum- 
ming probabilities for successive emissions using the Alterelli-Parisi splitting functions. Each 
term in the series that is resummed corresponds to an additional semi-classical radiation. The 
splitting functions only account for collinear emissions; to include soft emission, it is common 
either to impose some kind of angular ordering constraint to simulate soft coherence effects, or 
to use more sophisticated probability functions, such as Catani-Seymour dipoles [13]. Except 
for [H], none of these approaches has led to a resummation for event shapes beyond NLL. 

The approach to resummation of event shapes [T5] based on Soft-Collinear Effective The- 
ory (SCET) [16j [T7J [18] contrasts sharply with the semi-classical CTTW treatment. The 
most important conceptual difference is that effective field theory works with amplitudes, at 
the operator level, instead of probabilities at the level of a differential cross-section. Conse- 
quently, the resummation comes not from the exponentially decreasing probability for multiple 
emissions, but from a solution to renormalization group (RG) equations. 

The starting point for the effective field theory approach is the factorization formula for 
thrust in the 2-jet region, 

1 d(T 2 _ .A [a^A^AU 7^2 .A x/Q x C ( 1, .AXf~ Pl+Pr k 



(T dr 



H(Q\ fx) J dp 2 L dp 2 R dk J(pl fj) J{p% fx) S T (k } fi)S(r - - ^) , (4) 



where H(Q 2 , /i) is the hard function, J(p 2 , fx) the jet function, and Sr(k, fx) is the soft function 
for thrust. Q refers to the center-of-mass energy of the collision, fx is an arbitrary renormaliza- 
tion scale, and the born-level cross section do appears for normalization. A similar factorization 
formula was derived to study top quark jets in [19] . and then transformed into this form to 
study event shapes in [15] . Factorization properties of event shape variables were also studied 
in [201 [21]. The expression (jlj) is valid to all orders in perturbation theory up to terms which 
are power suppressed in the two-jet region r — > 0, 



da dcr 2 
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dr dr 

The key to the factorization theorem is that near maximum thrust, r reduces to the sum 
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of hemisphere masses 
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(6) 



where the two hemispheres are defined by the thrust axis n. Here, pf,(pfj) is the invariant 
mass of the energetic particles in the left (right) jet and kQ is the increase of the invariant 
mass on the two sides due to soft emissions. A more detailed interpretation of this formula 
can be found in ^ [191 E2] • 

The factorization theorem (jlj) makes it evident that the thrust distribution involves three 
different scales in the endpoint region. First of all, there are virtual effects arising in the 
production of the quark anti-quark pair at the hard scale /i^ ~ Q which are encoded in the 
hard function H(Q 2 ,fi). A second relevant scale is associated with the invariant mass of the 
two back-to-back jets, fi 2 ~ p 2 L + p\ ~ tQ 2 . In addition to these two external scales, a third, 
lower seesaw scale is encoded in the soft function fi s ~ k ~ tQ ~ fi 2 /fj>h- The effective theory 
treatment separates the effects associated with these three scales and makes transparent that a 
larger range of scales, and consequently a larger range of a s (ii) is being probed than is evident 
in either the fixed-order calculation or in the traditional NLL resummation. Large logarithms 
are avoided using RG evolution in the effective theory. Each of the three functions H, J and 
S is evaluated at its characteristic scale, and then evolved to a common scale /i. Solving the 
differential RG equations resums logarithms of the scale ratios. 

In Section El we provide the definitions of H, J and S in SCET. These functions can be 
calculated directly in SCET, or one can rewrite them in terms of matrix elements of QCD 
operators. For practical calculations, the definitions in QCD are often more suitable since the 
QCD Feynman rules are simpler. The hard and jet function appear in other processes and are 
known to two-loop order [231 EH EE]- Their RG-equations have been solved in closed form and 
the relevant anomalous dimensions are known at three- loop order [26]. With the hard and jet 
functions known, the only missing ingredient to resum the thrust distribution to next-to-next- 
to-next-to-leading logarithmic (N 3 LL) accuracy is the soft function S. Its one-loop expression 
was given in [T3] and in Section [3] we determine soft function to two loops. Its logarithmic part 
is obtained using RG invariance of the thrust distribution and the remaining constant piece 
by a numerical procedure. After plugging the solutions back into the factorization theorem 
(j4]), we obtain the result for the resummed distribution valid to N 3 LL. 

Next, we expand the effective theory result to fixed order in a s . The logarithmically 
enhanced terms which are determined in the effective theory dominate the thrust distribution. 
This is especially pronounced at NNLO: color structure by color structure, we find that the 
logarithmic terms are an excellent approximation of the full fixed-order result. This close 
agreement also provides an independent check on the NNLO calculation. After comparing the 
full fixed-order result to the logarithmic terms, we add the small difference between the two 
to our resummed result. By this matching procedure, we obtain a resummed result which is 
also correct to NNLO in fixed-order perturbation theory. 

In Section [H we fit the resummed matched calculation to aleph and opal data. We find 
a relatively small perturbative uncertainty on a s compared to previous event shape fits of 
the same lep data. In fact, the final statistical, systematic, perturbative and hadronization 
uncertainties end up being quite similar, all around 1%. At this point, we have the least 
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handle on hadronization effects, and these and other power corrections are explored in Section 
The conclusion contains a brief discussion of how the various uncertainties might be further 
reduced. 

2 Resummation of thrust in effective field theory 

The large logarithms in the thrust distribution dominate near the endpoint, r — * 0. This 
region of phase space corresponds to configurations with two back-to-back light jets. In this 
situation, the vector and axial-vector currents relevant to the production of the gg-pair are 
mapped onto the two-jet operators in SCET [27] 



where r = 7 M or T = 7^75 for vector or axial-vector currents respectively. Here, n is a light- 
like 4- vector aligned with the thrust axis and the composite fields Xn and Xn are the collinear 
quark fields in the n- and n-directions, multiplied by light-like Wilson lines [28] . The first step 
in the effective field theory calculation is matching to full QCD. This is done by calculating 
matrix elements in SCET and in QCD and adjusting the Wilson coefficients in the effective 
theory so that the matrix elements agree. Performing the matching on-shell, one finds that 
the relevant matching coefficient for the vector operator is given by the on-shell vector quark 
form factor. In a scheme with an anti-commuting 75, the Wilson coefficients of the vector 
and axial-vector operators are identical. Neglecting electro-weak corrections, the use of such 
a scheme is consistent in the endpoint region r — > 0. In this region, the two energetic quarks 
produced directly by the current always appear in the final state, so that the 75 matrices from 
the axial currents always appear in a single trace formed by the cut fermion loop. 

After normalizing to the tree-level cross section, the hard function H(Q, //) is given by the 
absolute value squared of the time-like on-shell form factor. Using the known two-loop result 
for the on-shell QCD form factor [291 USB EH [32], the hard function at two loops was derived 
in [23]. It satisfies the RG equation [26] 



2 = XfXx. 



(7) 



d 



H(Q 2 ,fi)= 2r cusp (a s )ln\ + 2 7 H (cO H(Q 2 ,fi) 




dln/i 



whose solution can be written as 



( 
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2A T {fi h ,n) 



H(Q 2 , fi) = H(Q 2 , fi h ) exp [4%h, n) - 2A H {^ h , a0] 



(9) 



Here, 
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and 
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The function Ar(v,fj) is defined as A H {y,n), but with ^ u replaced by r cusp . The solutions of 
the RG equations for the jet and soft function given below involve functions Aj(is, fi), A s {y, /i) 
which are obtained from (ITT]) by substituting 7j,7s for 7^. It is straightforward to expand 
S(v, /i) and v4#(i/, /1) perturbatively in a s (z/) and a s {ii) given the expansions of r cusp (o;) and 
7 (a). The explicit expansions can be found in [24J. 

In SCET the jet function is given by the imaginary part of the collinear quark propagator, 



[n ■ p) 7r 



rf 4 xe-^(0|T^n(x)|xn(0)MO) 



5{p 2 ) + 0{a s ) (12) 



and thus vanishes for p 2 < 0. The jet function was calculated at one loop in [33] and at two 
loops in [23j. To evaluate the function perturbatively, it is convenient to rewrite the collinear 
quark propagator in terms of QCD fields. One finds that the jet function is obtained from 
the quark propagator in light-cone gauge. The jet function satisfies a RG equation which is 
non-local in p 2 |23j, 

MPV) [„ ,...P> „. l J(pV) + 2r „^ V<^>- J <^> . (13) 

Jo P - Q 



dln/i 



-2r cusp In — - 27j 
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From the divergent part of the form factor at three loops [32j and the NNLO Altarelli-Parisi 
splitting functions [31] the jet anomalous dimension 77 was derived at three loops in [21] and 
is given in Appendix |A] Although the RG equation is non-local in p 2 , it is local in \i and can 
be solved using Laplace transform techniques. The result is 



~ 1 fp 2 \ Vi e~ jEVj 

J(p ,») = exp[-AS(fi j ,n) + 2Aj(fi j ,fi)]j(dr h ,fj. j )-^ , (14) 

where rjj = 2A-p(fij, fi). The function j(L,fi) is the Laplace transform of the jet function. Its 
definition and explicit form are given in Appendix [Bl To any given order in perturbation 
theory, j(L, /x) is a polynomial in the variable L so that the derivatives with respect to rjj in 
( |T4l) can be performed explicitly. 

The thrust soft function is defined as a matrix element of Wilson lines along the directions 
of the energetic quarks, 



S T (k) = I (X\Y^Y n \0)\ 2 6(k - n -p Xn -n- PXn ), (15) 

where 



x 



Y n = P exp ^ig J dtn ■ A s (tn)^j . 



(16) 



This Wilson line describes the Eikonal interactions of soft gluons with the fast moving quark, 
and P x n {pxn) ls the sum of the momenta of the soft particles in the n- hemisphere (n- 
hemisphere). The variable k measures the change in the invariant mass due to soft emissions 
from the two jets. At the leading power, the mass in the n-hemisphere is given by 

M 2 n = (p n + pxj ^p 2 n + Q(n- p Xn ) , (17) 
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where p n denotes the total collinear momentum in the hemisphere. Note that the soft function 
vanishes for negative argument. 

Like the jet function, the soft function can be calculated order-by-order in perturbation 
theory. The one- loop soft function was derived in [15] from results of [35] ; it was also calculated 
directly in SCET [22]. The two- loop soft function will be determined below. 

The factorization theorem (T4J) and the fact that the thrust distribution is independent of 
the renormalization scale fi implies that the soft function fulfills the RG equation 



4r 



dS T (k,n) 
din /i 

and that, to all orders 



k 



cusp(as) In 27 (a s 



S T (k, 11) — 4T 



cusp \C^s j 







k ^ S T (k^)-S T (k'^) 



7 s = 1 H -2j J . (19) 

This relation was checked to one loop in [15] (with a different convention for 7#), and here we 
use it to determine the two- and three- loop soft anomalous dimensions. Similar to (|14p . the 
solution for the soft function is 



S T {k, n) = exp [4S (/j, s , ji) + 2A s (fi s , //)] sr(^ s 



1 fk \ Vs e~ lEVs 



(20) 



with r] s = — 4Ar(/i s , /i). From the linearity of As in 75 it also follows that As = Ah — 2Aj. 

The convolution integrals in (J4J) can be done analytically once the solutions (fl4l) and ([20 
are put back into the factorization theorem. The thrust distribution becomes 



1 da 2 



x 



do dr 

Q2 \ -2j4 r (Mh,M) 



exp [4S(ii h ,ii) - 2A H (n h ,n) - 8S(fij,fi) + 4Aj(fij,fi) + 45(/x s ,/i) + 2A s (fj, s , fx)] 



i 2 



H(Q ,fi h ) j(d2 Vj ,Vj) s T (d VB , fi 8 ) 



1 /rQ 2 



2 \ 2% 



rQY h e 



(21) 



Using the relations, 



Ar(vi, 1^2) + A r (fi2, fa) = A T (fa, fa) 



fa 



5(//i,// 2 ) + S(u, 2 ,fa) = S(fa,fa) +ln— A T {^ 2 ,fa) 

fa 



and 



f{d v )X* = XVf(\nX + d v ) 



the expression (j2"Tj) simplifies to 



(22) 



(23) 



— -p 1 - = exp [4S(fi h ,fij) + 4S(u, s , fij) - 2A H (fa,fa) + AAj(fij, fi s )\ ( % 
cr ar \n h 



2\ -2A r (w>,Mj) 
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order 


■p 

1 cusp 


ryH/J/S 


tl, J, s T 


a 
P 


fixed-order 
matching 


logarithmic 

accuracy 


l st order 


2-loop 


1-loop 


tree 


2-loop 




NLL 


2 nd order 


3-loop 


2-loop 


1-loop 


3-loop 


LO 


NNLL 


3 rd order 


4-loop 


3-loop 


2-loop 


4-loop 


NLO 


N 3 LL 


4 th order 


4-loop 


3-loop 


3-loop 


4-loop 


NNLO 


N 3 LL 



Table 1: Definition of orders in perturbation theory 



xH(Q 2 ,fi h ) 



■I 1 , r. 



1 2 



-~fEV 



T( V ) 



(24) 



with r] = 4A-p(nj, /j, s ). From this final result we can read off the canonical relations among the 
hard, jet, and soft matching scales and the physical scales Q and p ~ y/rQ: 



Hh = Q, Hj = VtQ , Hs = tQ 



(25) 



Note that the arbitrary reference scale \x has dropped out completely. 

For the a s fits, we need the differential thrust distribution integrated over each bin. The 
integral of the thrust distribution can be evaluated analytically, since the derivatives with 
respect to rj in ( l24|) commute with the integration over r. The resulting expression is 



R?.(t) 



X 



T 1 da 2 



-2A F (ii h ,tJ,j) 



T dr' = exp [4S(ji h ,Hj) + 4S(/i s ,Hj) - 2A H (p, h ,n s ) + 4A J (// i , /i s )] 



H(Q 2 ,fi h ) 



A 4 ? 



-i 2 



tQ 

v-s) r(77 + 1) 



(26) 



Note that the integral is performed for fixed \ij and /x s , that is, before setting them to their 
canonical r-dependent values. In this way, large logarithms are removed in the observable of 
interest, not for some intermediate expression. 

Different definitions of logarithmic accuracy are commonly used in the literature. Before 
proceeding further, we now show which logarithms are included at a given order in our calcu- 
lation. We use renormalizat ion-group improved perturbation theory, in which logarithms of 
scales are eliminated in favor of coupling constants at different scales which are counted as 
small parameters of the same order 



In 



da 



2tt 



1 



,{u) A) \a,(p) a„(u) 

The expansion of the Sudakov exponent ( fTOl) then takes the form 

1 



+ 



(27) 



3(v, fi) 



/i(r) + / 2 (r) + a s {v)h{r) + a s (v) 2 f 4 (r) + ... 



(28) 
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where r = a s (fi)/a s (u). The explicit expressions for the functions f\ to / 4 needed for our 
calculation are given in [24J. The leading-order term a° s in renormalization group improved 
perturbation theory involves the functions fi and f 2 , which depend on the one and two-loop 
cusp anomalous dimension. To make contact with the literature, we can expand oi s (n) around 
fixed coupling a s = a s (z/). The result takes the form 

S(v,[i) = Lg x (a s L) + g 2 (a s L) + a s g 3 (a s L) + a 2 s g A (a s L) + ... , (29) 

with L = ln(/x/z/). LL resummations include only g±, NLL also g 2 and so forth. When 
rewriting (1281) in the form (1291) . the expansion of fi contributes to the functions gj with 
j > i so that there is a one-to-one correspondence between the order in renormalization 
group improved perturbation theory and the standard logarithmic accuracy. Note that the 
higher order terms to (128]) and ff29l) are suppressed by explicit factors of a s . The missing 
pieces in the integral R 2 {t) at N 3 LL are suppressed by a 3 so that the missing logarithms are 
a 3 x a n ln 2n r = a k ln 2fc_6 r for the default scale choice. In particular, at order a 3 the N 3 LL 
result includes everything except for the constant term in R 2 {t) which does not contribute to 
the thrust distribution. 

In Table [U we list the ingredients to obtain (1261) to a given accuracy. The necessary anoma- 
lous dimensions and the results for the functions H, j and St are provided in Appendix lAl 
Everything in the table except for the four-loop cusp anomalous dimension and the constant 
part of the two-loop soft function are known. We estimate the former using the Pade approxi- 
mation T 4 = Tl/T 2 [36] and determine the latter numerically in the next section. Rather than 
specifying both the accuracy of the resummation and the order to which we match to the fixed 
order result, will will in the following simply refer to the definitions of 1 st , 2 nd , 3 rd and 4 th order 
as given in Table [U Note that the difference between 3 rd and 4 th order, as we have defined 
them, is only the inclusion of NNLO matching corrections, but the logarithmic accuracy stays 
the same. 

3 Resummation vs. fixed order 

In this section, we compare the resummed expression, valid in the endpoint region r — > to 
the fixed-order expression, which is valid away from the endpoint. The resummed expression, 
when expanded to fixed order, must reproduce the r = singularities of the fixed-order 
calculation. This observation can be used to extract numerically the constant part of the 
two-loop soft function. Then, by including the difference between the expanded resummed 
expression and the fixed-order expression, we derive the final matched distribution. 

The fixed-order thrust distribution has been calculated to leading order analytically and 
to NLO and NNLO numerically. For the scale choice fi = Q, the result is usually written in 
the form 

^^m A ^(^BiT )+ (^C {T) + ---, (30, 

where we have suppressed the argument of the coupling constant, a s = a s {Q). Throughout 
the following analysis, we use an analytical form for A(t), a numerical calculation of B{r) 
using the program event2 [37J with 10 10 events and a numerical calculation of C(r) that was 
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Figure 1: A comparison of the full fixed-order calculations and the fixed-order expansion of 
the resummed distributions from the effective field theory. The light-red areas in the NNLO 
histogram are an estimate of the statistical uncertainty. 



generously provided by the authors of [2J. A value of yo = 10~ 5 for the infrared cut-off was 
used in the calculation of the NNLO histograms, see |38j. 

The resummed differential thrust distribution in the effective theory is given in Eq. (|24p . 
To compare with fixed-order results ( l30l) . we set all scales equal Hh = jij = fi s = Q. Doing so 
switches off the resummation: all evolution factors, such as S(fih,fij) and Au^hi vanish 
in the limit of equal scales. Before taking the limit 77 = 4Ar(fij, /i s ) — > 0, we expand the kernel 
in (12^1) using the relation 



1 1 n n 

-r- = V) + E^r 

1 n=0 



ln n r 



(31) 



and evaluate the derivatives with respect to 77 using the explicit expressions for j and s. The 
result is a sum of distributions 



a dr S ^ )D ' + (I) + (s) + (I) + " " (32) 

The coefficients D$, Da Db and Dc are given in Appendix O Away from r = 0, the 8- 
function terms can be dropped and the plus-distributions reduce to their argument functions, 
[Dx(r)] + = Dx{t). Since the effective field theory resums the large logarithms of the fixed- 
order distribution, there should not be any 1/r singularities in A, B, or C which are not 
reproduced in Da, Db and Dc respectively. This was shown analytically for the A function 
in [15] . It is demonstrated numerically for A, B and C in Figured! In fact, the figure shows 
that even at moderate r, the thrust distribution is dominated by the singular terms. Note that 
the lowest three bins of the numerical result for C are above the effective theory prediction. 
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This is due to numerical difficulties in the fixed-order code used to evaluate C and will be 
explored in more detail below. 

The SCET expression (I24p for the thrust distribution is valid as r —>■ 0, that is, in the 
2-jet region. One could perform resummation also for terms which are power suppressed 
in this limit, by including operators with additional fields or derivatives into the effective 
theory [39], HQ]. However, since these terms are power suppressed it is sufficient to include 
them at fixed order. To do so, we simply subtract the singular terms from the fixed-order 
expression. The remainder is 



r(r) 



1 f da dcr 2 



o"o Vdr dr 



(£) [A{T) ~ DA{T)] + (i) 2 [B{T) ~ Db{t)] + (i) 3 [C{T) ~ Dc{T)] + 



Including the matching contribution, the thrust distribution becomes 

1 da 1 dcr 2 
a dr 



+ r r 



(33) 



(34) 



o- dr 

With the inclusion of r(r), our result not only resums the thrust distribution to N 3 LL, but is 
is also correct to NNLO in fixed-order perturbation theory. 

Now let us turn to the two-loop soft function. Its RG equation together with the anomalous 
dimensions determine the logarithmic part of s^, but the constant part 

2 



sr(0,/z) 



°f clc F + C A c s 2 Ca + T F n f c: 



'2,n, 



(35) 



cannot be obtained in this way. We will determine the constant from the requirement that 
the integral over the thrust distribution reproduces the total hadronic cross section 



°had 



/123 \ 3 

C F C A ( — - 44C 3 J + C F T F n } (-22 + 16Ca) - C 2 F - 



■ (36) 



10 



100 F 



50 



-50 



-100 I- 
0.000 



0.001 



0.002 



0.003 



0.004 



0.005 



l-T 



Figure 3: Extraction of the two-loop constants in the soft function. The points correspond to 
the value of an infrared cutoff applied to the fixed-order calculation. The lines are interpo- 
lations among the points from r = 0.001 to r = 0.003 extrapolated to r = to extract the 
constants. From top to bottom, the curves are the Cp, Ca and n/ color factors. 



Plugging (PJ) and (J33D into fl34]) we find 



Chad 



D 5 + 



dr r(r) 



i + — + M 

7T V47T 



317.5 + cf + 4 / [B(t) - D b {t)} dr \ , (37) 



where 317.5 comes from setting rif = 5 in D$ (see Appendix C). Since we know separately 
the color structures for B (numerically) and D% (analytically), as shown in Figure [2j we 
can perform this integral numerically and then extract cf by comparing to ( 1361) . Although 
the difference B{r) — Db(t) is integrable as r — > both of these functions are separately 
divergent. To have numerically stable results, we impose an infrared cutoff To on the integral 
and interpolate to r = 0. We do this in discrete steps by dropping the lowest bins in the 
B(t) distribution which was generated with the event2 program. The convergence and 
interpolation are shown in Figure [3J We find 



"2 CF 



58 ±2. 



-60 + 1 



43 + 1 



(3? 



These constants were explored previously in [12]. Lacking the form of the divergences near 
r = 0, these authors had to fit for the shape of the curve as well as the constants, leading to 
results with much poorer accuracy. A comparison with the results of rj2] is given Appendix 


We now have all the necessary perturbative input at hand to evaluate the thrust distribu- 
tion and to extract a s . Before doing so, we compare the recent NNLO fixed-order results in 
detail to the singular terms predicted using the effective theory. In Figure H] the contribution 
of the six color structures which appear at to C{r) and Dq{t) are plotted. The color 
structure of the NNLO coefficient C has the form C = C F (N 2 d + C 2 + l/iV 2 C 3 + Nn f C A + 
rif/NCs + njCe) and the plot shows the six parts, with the prefactors evaluated for iV = 3 
colors and nj = 5 quark flavors. The figure shows that the singular terms (blue lines) are a 
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Figure 4: Contributions of different color structures to the three-loop coefficients of the thrust 
distribution. The plots show a comparison of our result for the singular terms encoded in Dq 
(blue lines) with the numerical evaluation of the full coefficient C (red histograms) [38J. The 
light-red areas are an estimate of the statistical uncertainty. 

good approximation to the full result (red lines) for each color structure. What is surprising is 
that they seem to agree well almost everywhere. One consequence of this is that the matching 
to the NNLO fixed-order distributions will have a small effect. The dominance of the loga- 
rithmically enhanced terms, even at moderate r, strongly suggests that resummation would 
indeed lead to a significant improvement in perturbative accuracy. The close agreement also 
provides a verification of the fixed-order result. Because the same numerical code is used for 
many other NNLO observables, such an independent check is certainly welcome. 

As we observed earlier, the lowest three bins of the NNLO fixed-order result of [38] are 
higher than the singular terms obtained with the effective theory, see Figure [U The excess 
at small r seen in Figure [1] is barely noticeable in Figure HI because we have multiplied the 
distributions by r which de-emphasizes the small-r region. To analyze this region in detail, 
we plot the distribution as a function of In r in Figure [5j For very small r, the full result 
should reduce to the singular terms derived in the effective theory. However, this region is 
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Figure 5: Contributions to the three-loop coefficients of the thrust distribution. The plots show 
a comparison of our result for the singular terms (blue lines) with the numerical evaluation 
of the full result (red histograms) [38]. The dotted, dashed and solid lines correspond to an 
infrared cut-off yo = 10 _5 ,10~ 6 and 10 -7 , see |38j . The light-red areas are an estimate of the 
statistical uncertainty. 



very challenging for the numerical integration. The numerical results are shown in red in 
Figure and the light-red bands are the statistical uncertainty from the numerical NNLO 
calculation. The three red lines correspond to different values of an infrared cutoff, which 
is imposed when generating events [3B|. The agreement is good, except for the two leading 
color structures. The authors of [38] are aware of the problem [H]. For the extraction of a s , 
the region of very small r will not be used, so these numerical difficulties are not critical for 
present purposes. 
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Figure 6: Convergence of resummed and fixed-order distributions, aleph data (red) and opal 
data (blue) at 91.2 GeV are included for reference. All plots have a s {mz) = 0.1168. 



4 a s extraction and error analysis 

In this section we now use our result for the thrust distribution to determine a s , using lep data 
from ALEPH [42 J and OPAL |43J. Before performing the fit, let us compare the perturbative 
expansion with and without resummation. The result at Q = 91.2 GeV is shown in FigureE] 
side-by-side with the fixed-order expression. We use the same value a s (mz) = 0.1168 for both 
plots and have set the scales /i^, /i-,- and fi s to their canonical values ff25l) . For reference, we 
also show the aleph and opal data. The curves for the fixed-order calculation correspond to 
the standard LO, NLO, NNLO series; for the effective field theory calculation, the orders are 
defined in Table [TJ It is quite striking how much faster the resummed distribution converges. 
In fact, it is hard to even distinguish the higher order curves after resummation, except in 
the region of very small r, where the distribution peaks. The peak region is affected by non- 
perturbative effects, as will be discussed in the next section, but it will not be used in the 
extraction of a s . The region relevant for the a s extraction is shown in the lower two plots. 
The value of a s (mz) = 0.1168 we use in the plots corresponds to the best fit value in the range 
0.1 < t < 0.24 for the aleph data set. However, the plot makes it evident that the extracted 
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Figure 7: Relative error for best fit to aleph data at 91.2 GeV. The inner green band includes 
only statistical uncertainty, while the outer yellow band includes statistical, systematic and 
hadronization uncertainties. The solid line is fit to 0.1 < 1 — T < 0.24 giving a s (mz) = 0.1168 
while the dashed line is fit from 0.08 < 1 — T < 0.3 giving a s (mz) = 0.1171. The smaller fit 
range is used for the error analysis because it has been previously studied in [5]. 



a s value will not change much beyond first order. A fit to the NNLO fixed-order prediction 
gives a s {m z ) = 0.1275. 

The aleph and opal collaborations have published analyses of the lep 1 and higher 
energy lep 2 thrust distributions. To fit a s we calculate the thrust distribution integrated 
over each bin measured in the experiments. The resummed contribution in a given bin is 
obtained as ^(tr) — RiijL) using Eq. fl26|) for the bin with n < r < tr. For the matching 
contribution, we integrate analytically the Da(t), Db(t) and Dq{t) functions and subtract 
them from the analytic integral of A(t) and the appropriately binned numerical distributions 
B(t) and C(r). 

A problem we encounter when trying to extract a s is that the experiments have published 
statistical, systematic, and hadronization uncertainties for each bin, but have not made the 
bin-by-bin correlations public. Without this information, we proceed with a conservative 
approach to error estimates: to extract the default value of a s , we perform a x 2_ fit to the 
data including only statistical uncertainties. We then use the systematic and hadronization 
errors on a s obtained in previous fits to aleph [5] and opal [13] data. In these papers fits 
to a s were performed which included the correlation information. To be able to use their 
values, we perform our fits using exactly the same fit ranges as used in these papers. This 
is not entirely optimal, since the experimental systematic error will depend somewhat on the 
theoretical model used in the fit. Our resummed calculation is valid in a wider range of r than 
the predictions used in [51 03], so one could use data closer to the peak, where the statistics 
are higher and resummation is more important. In a future analysis, the fit range could be 
optimized to minimize the total error after folding in the proper correlations. 

In Figure El we plot the relative statistical and total experimental uncertainty as a function 
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Figure 8: Perturbative uncertainty at Q = 91.2 GeV. The first four panels show the variation 
of the matching scale, the hard scale, the jet scale, and the soft scale. Each of the scales is 
varied separately by a factor of two around the default value. The last two panels show the 
effect of simultaneously varying the jet- and soft scales, see text. The lep 1 aleph data is 
included for reference. All plots have a s (mz) = 0.1168. 

of t and compare to the best fit result. We find that the extracted value is fairly insensitive to 
the fit range. In fact, going from the standard range (solid line) to the larger region (dashed 
lines) changes the best-fit value of a s {mz) by less than 0.3%, from 0.1168 to 0.1171. 

Next, we consider the perturbative theoretical uncertainty. In the effective field theory 
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analysis, four scales appear: the hard scale /i/j ~ Q, the jet scale /ij ~ \/tQ, the soft scale 
fi s ~ tQ, and the scale fi m at which the matching corrections are added. In the matching 
corrections the physics associated with the hard, jet and soft scales has not been factorized, so 
it is not obvious which value of /i m should be chosen. We follow standard fixed-order practice 
and choose /i m = Q as the default value. Our result is independent of these scales to the order 
of the calculation: the change in the result due to scale variation can thus be used to estimate 
the size of unknown higher order terms, of 0(a^) for our final result. 

We show the results of varying each of the four scales up and down by a factor of 2 in the 
first four panels of Figure [HI The results converge nicely, with the dominant uncertainty coming 
from the soft scale variation. This is expected, as the soft scale probes the lowest energies 
and therefore the largest values of a s . In fact, it is a critical advantage of the effective theory 
that the soft scale can be probed explicitly - the fixed-order calculation has access to only one 
scale and assuming fi ~ Q may therefore underestimate the perturbative uncertainty. From 
the first panel in Figure [HJ it is clear that the extraction of a s is almost completely insensitive 
to the scale at which the fixed order calculation comes in. Again, this is in contrast to a pure 
fixed-order result. The matching scale variation is so small because the matching correction 
itself is small, as we saw in Figures [fl [2J HJ and0 

Figure [HJ shows the effect of varying the jet and soft scales separately by factors of two: 
\\frQ < fij < 2y/rQ and \tQ < jj, s < 2rQ. While a factor of two may seem reasonable for a 
fixed order calculation (although as we have already observed, the thrust distribution probes 
scales tQ <C Q), from the effective field theory point of view it makes little sense to vary the 
soft and jet scales separately. In doing so, one can easily have fij < fi s or fih < fij which is 
completely unphysical. Instead, for the error analysis we will use two coordinated variations. 
First, a correlated variation holding fij/fi s fixed: 

\Lj — > cy/rQ, jj, s — > ctQ, i < c < 2 . (39) 

This probes the upper and lower limits on fij and fi s , but avoids the unphysical region. Second, 
an anti-correlated variation, holding fi 2 /(Qfi s ) fixed: 

fi 2 — > aQ 2 r fi s — > aQr, < a < V2 . (40) 

v2 

This is independent from the correlated mode but again avoids having fij < ft s . The uncer- 
tainty resulting from these two variations is shown in the last two panels of Figure [HJ 

To estimate the total perturbative uncertainty on the extracted value of a s , we use the 
uncertainty band technique proposed in jH] and adopted both by aleph [32] and opal [32] 
as well as in the recent fit of NNLO results to aleph data [5]. The result is shown in Figure [9j 
In short, the theoretical uncertainty is determined as follows: one first calculates a s (mz) using 
a least-squares fit to the data with all scales at their canonical values and without including 
any theoretical uncertainty in the x 2- f unc tion. Then each scale is varied separately, holding 
ots{^z) fixed to its best-fit value. These produce the curves in FigureEl Next, the uncertainty 
band, the yellow region in Figure El is defined as the envelope of all these variations. Finally, 
the scales are returned to their canonical values, and the maximal and minimal values of a s are 
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Figure 9: Uncertainty bands for various scale variations. The band in the first panel is deter- 
mined entirely by scale variations. The second panel shows an alternative way of estimating the 
perturbative uncertainty using an educated guess of the uncalculated higher order coefficients, 
as described in the text. 
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Q 


91.2 


133 


161 


172 


183 


189 


200 


206 


AVG 


nt range 


0.1 
0.24 


0.06 
0.25 


0.06 
0.25 


0.04 
0.2 


0.06 
0.25 


0.04 
0.2 


0.04 
0.2 


0.04 
0.2 




xVd.o.f. 


32.5/13 


7.7/4 


3.3/4 


10.3/4 


3.6/4 


0.9/4 


24.6/4 


4.0/4 


- 


stat. err. 


0.0001 


0.0037 


0.0070 


0.0080 


0.0043 


0.0022 


0.0023 


0.0023 


0.0010 


syst. err. 


0.0008 


0.0010 


0.0010 


0.0010 


0.0011 


0.0010 


0.0010 


0.0010 


0.0010 


hadr. err. 


0.0019 


0.0014 


0.0012 


0.0012 


0.0011 


0.0011 


0.0010 


0.0010 


0.0012 


pert. err. 


+0.0013 
-0.0017 


+0.0012 
-0.0016 


+0.0015 
-0.0020 


+0.0006 
-0.0009 


+0.0010 
-0.0013 


+0.0011 
-0.0015 


+0.0010 
-0.0014 


+0.0009 
-0.0012 


0.0012 


tot. err. 


0.0026 


0.0043 


0.0074 


0.0082 


0.0047 


0.0030 


0.0030 


0.0029 


0.0022 


(Pade x 2) 


0.0004 


0.0004 


0.0003 


0.0004 


0.0004 


0.0004 


0.0004 


0.0003 




a s (m z ) 


0.1168 


0.1183 


0.1263 


0.1059 


0.1160 


0.1203 


0.1175 


0.1140 


0.1168 


PYTHIA 


0.1152 


0.1164 


0.1248 


0.1028 


0.1146 


0.1177 


0.1151 


0.1119 


0.1146 


ARIADNE 


0.1169 


0.1181 


0.1264 


0.1047 


0.1164 


0.1197 


0.1170 


0.1135 


0.1164 



Table 2: Best fit to aleph data. The row labeled (Pade x 2) is an alternative measure of 
perturbative uncertainty as described in the text. It is not combined into the total error. The 
rows labeled pythia and ARIADNE give the value of a s after correcting for hadronization and 
quark masses using pythia or ARIADNE. 



determined which allow the prediction to remain within the uncertainty band. An important 
feature of this approach is that the data enters only in the determination of the best fit a s and 
the fit region; the perturbative uncertainty is determined purely from within the theoretical 
calculation. Separating the theoretical and experimental errors in this way makes it much 
easier to average a s results obtained from different data sets, since they suffer from the same 
theoretical uncertainty. 

The purpose of scale variations is to estimate the effect that a higher order perturbative 
calculation would have on a distribution. This is justified by arguing that any scale variation 
can be compensated by terms at one order higher in a s , thus it should give a reasonable 
estimate of these higher order terms. However, as we have seen, the amount by which we vary 
the scales is arbitrary, and the traditional factor of 2 in the variation is both problematic for 
the jet and soft scales and seems to overestimate the uncertainty. The distribution determined 
by the effective theory at one higher order depends on only a handful of numbers: the beta- 
function coefficient /?4, the anomalous dimensions T^^^j^ and the constants in the hard, 
jet and soft functions, c^,c^,c^, see Appendix IA1 for the subscript conventions. Thus in the 
effective field theory, there is a straightforward way to estimate the effect of higher orders: 

one simply varies these coefficients. For example, we can estimate their size using a Pade 

r 2 

approximation: T n +i = ic^- 2 -. This should reasonably span likely values for what a higher 

J- 71 — 1 

order perturbative calculation would provide. We show the variations corresponding to c = 2 
and c = 5 in the bottom panel of Figure M, which are labeled Pade x 2 and Pade x 5 
respectively. In each case we scan over the signs for the various coefficients to find the largest 
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Q 


91.2 


133 


177 


197 


AVG 


fit range 


0.05-0.3 


0.05-0.3 


0.05-0.3 


0.05-0.3 




X 2 /d.o.f. 
stat. err. 


149.9/5 
0.0001 


17.0/5 
0.0038 


1.7/5 
0.0033 


18.3/5 
0.0014 


0.0014 


syst. err. 


0.0011 


0.0054 


0.0028 


0.0013 


0.0013 


hadr. err. 
pert. err. 
tot. err. 


0.0031 

+0.0014 
-0.0018 

0.0037 


0.0024 

+0.0011 
-0.0015 

0.0072 


0.0021 

+0.0009 
-0.0013 

0.0049 


0.0019 

+0.0011 
-0.0014 

0.0030 


0.0019 
0.0013 
0.0030 


(Pade x 2) 


0.0004 


0.0003 


0.0003 


0.0003 




a s (m z ) 


0.1189 


0.1165 


0.1153 


0.1189 


0.1189 


PYTHIA 


0.1143 


0.1142 


0.1134 


0.1173 


0.1173 


ARIADNE 


0.1163 


0.1160 


0.1151 


0.1189 


0.1189 



Table 3: Best fit to opal data. 



variations. Even for c = 2, the fifth order coefficients come out quite large, for example, 
T.4 pa ±2 x 10 4 and 04 ~ ±3 x 10 5 . Nevertheless, the uncertainty is still significantly smaller 
than what we found using scale variation. We find that Pade x 2 gives 5a s (mz) ~ 0.0003 in 
contrast to errors around 5a s (mz) ~ 0.0012 from the scale variations. Although the higher 
order constants are unknown, one might try to estimate them in more sophisticated ways, for 
example, by computing the dominant diagrams. In the end, we will not use this new method 
for the final error estimates, but we present the resulting uncertainties in Tables [2] and [3] 
for completeness. They include a scale variation in the matching correction because this is 
independent of the resummed distribution. 

For each of the energies in the aleph and OPAL data sets, we perform a least-squares fit 
using the experimental statistical errors. The statistical uncertainty on a s is calculated from 
the variation in y 2 , the perturbative uncertainty is calculated using the uncertainty band 
(with scale variations), and systematic uncertainties are taken from [5] and [13], as discussed 
above. We include the non-perturbative hadronization uncertainties from these papers, but 
do not include the corresponding hadronization corrections. We will discuss hadronization 
and other power-suppressed effects in detail in the next section. The fit results are given in 
Tables M and GJ 

To combine the results from different energies, we compute a weighted average, a s = 
WiOtf^ . The weights Wi are determined by minimizing the uncertainty a 2 = . Wi Wj cov(i, j) 
Given that we don't know the exact correlations, we set 

cov(z, j) = 2 5 id + og, og> + + °SW2t ■ (41) 

That is, we assume uncorrelated statistical errors and 100% correlation for the systematic, 
hadronic and perturbative uncertainties at different energies. Because the correlated uncer- 
tainties are dominant, naively minimizing the uncertainty can in some cases be lead to negative 
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Figure 10: Best fit values for a s (mz)- From right to left the lines are the total error bars at 
each energy for first order, second order, third order and fourth order, as defined in the text. 
The bands are weighted averages with errors combined from all energies. 



weights. This happens when combining the OPAL results in the above way. We eliminate 
these solutions by imposing Wi > 0, after which the best value from OPAL is obtained by 
assigning 100% weight to the highest energy measurement which has the smallest systematic 
uncertainty. The result obtained after combining aleph and opal results individually is given 
in the last column in Tables [2] and [31 Finally, we combine the aleph and opal results to an 
overall average. In this case, we assume that the systematic uncertainties are completely corre- 
lated between the individual energy results from each experiment, but neglect the correlations 
between the systematical uncertainties among the two experiments. For the hadronization 
and perturbative error, we assume 100% correlation. Proceeding in this way, we find 

a s (m z ) = 0.1172 ± O.OOlO(stat) ± 0.0008(sys) ± 0.0012(had) ± 0.0012(pert) 

= 0.1172 ±0.0022. (42) 

This result is close to the PDG world average a s (mz) = 0.1176 ± 0.0020 and has similar 
uncertainties. 

Our calculation does not include hadronization corrections and neglects quark masses. If 
we estimate their effect using pythia, the central value shifts to a s (mz) = 0.1150, while 
correcting with ARIADNE gives a s (mz) = 0.1168. We observe that the difference we find 
between pythia and ARIADNE is larger than the hadronization uncertainty in our average, 
which is based on aleph and opal studies. Correcting our higher order perturbative result 
with a tuned leading-order Monte Carlo shower is problematic, so this difference should be 
interpreted with caution. Various issues associated with hadronization corrections will be 
discussed in detail in the next section. 

It is interesting to repeat the fit order by order. This is done in Table H] and displayed 
graphically in Figure [TDJ The figure shows that the results found at different energies are 
consistent and illustrates the reduction of the uncertainty when including higher order terms. 
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ALEPH 




LEP 1 +LEP 2 


LEP 1 


order 


a s 


total err 


pert, err 


a s 


tot. err 


pert, err 


l st order 


0.1142 


0.0297 


0.0296 


0.1142 


0.0297 


0.0296 


2 nd order 


0.1152 


0.0068 


0.0064 


0.1166 


0.0071 


0.0068 


3 rd order 


0.1164 


0.0033 


0.0027 


0.1166 


0.0037 


0.0031 


4 th order 


0.1168 


0.0022 


0.0012 


0.1168 


0.0026 


0.0015 



OPAL 




LEP 1 + LEP 2 


LEP 1 


order 


a s 


total err 


pert, err 


a s 


tot. err 


pert, err 


l st order 


0.1190 


0.0305 


0.0304 


0.1190 


0.0305 


0.0304 


2 nd order 


0.1198 


0.0076 


0.0070 


0.1205 


0.0081 


0.0074 


3 rd order 


0.1194 


0.0040 


0.0029 


0.1194 


0.0047 


0.0034 


4 th order 


0.1189 


0.0030 


0.0013 


0.1189 


0.0037 


0.0016 



Table 4: Best fit values and uncertainties at different orders, as defined in Table [H 

5 Non-perturbative effects and quark mass corrections 



Let us now turn to two power suppressed effects which we have so far neglected in our analysis. 
The first is hadronization: the effective theory calculation corresponds to a parton-level distri- 
bution, while the experiment measures hadrons. Secondly, we have neglected quark masses in 
our calculation. Because thrust is an infrared-safe observable, both corrections are expected 
to be small, however they may not be negligible. 

Most of the previous determinations of a s have used Monte Carlo event generators to 
correct the parton-level predictions for hadronization effects and estimate the hadronic un- 
certainty by comparing the output of different generators. In particular, aleph [42] . opal 
[45] and the recent NNLO analysis [5] all use pythia to obtain their default hadronization 
corrections and then compare to herwig and ARIADNE to obtain the associated uncertainty. 
It turns out that the largest differences generally occur between pythia and ARIADNE |43j . 
even though ARIADNE uses pythia to calculate hadronization. 

We include in Tables [2] and [3] the best fit values of a s obtained after correcting the data bin- 
by-bin for hadronization and b- and c-quark masses using pythia v. 6. 409 [35], with default 
parameters, and with ARIADNE v. 4. 12 [36], using the aleph tune. Correcting with ARIADNE 
has quite a small effect on the values of a s . Moreover, ARIADNE always gives a larger value 
of a s than pythia. If the central values are taken after the ARIADNE corrections, they agree 
quite closely with a fit to the parton level distributions, that is, without any hadronization. In 
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addition, we also used the new sherpa dipole shower [47] for hadronization and find results 
similar to ARIADNE. 

Relying on the Monte-Carlo generators for hadronization is clearly not ideal, since they 
have been tuned to the same lep data we are trying to reproduce! The situation is especially 
problematic when trying to correct our resummed distribution. The Monte Carlo generators 
are all based on the parton-shower approximation, which only sums the leading Sudakov 
double logarithms and part of the next-to- leading logarithms. In contrast, our distribution is 
correct to N 3 LL and to NNLO in fixed-order perturbation theory. By tuning to data, part of 
the missing higher order perturbative corrections get absorbed into the hadronization model. 
An obvious way to avoid this problem would be to include the higher order corrections into 
the Monte Carlo codes, but needless to say, no such generator yet exists (although, see [39], [40] 
for an approach to improving generators based on the same effective field theory ideas we are 
using) . 

As shown in Figure [12] pythia agrees with the aleph data better than our 4th order 
resummed and matched theoretical calculation. How is this possible in a leading-log shower 
with leading-order matrix elements? The answer is that part of what is being tuned to data in 
the Monte Carlo program is not just the hadronization model but also some kind of unfaithful 
imitation of subleading-log resummation. This is demonstrated in Figure d2] where pythia 
is run at the parton and hadron level and compared to the 1st order and 4th order resummed 
matched distributions in the effective field theory. Even at the parton level, pythia agrees 
more with the 4th order than the 1st order. Moreover, the hadronization corrections provide 
something like a shift in the distribution, but cannot explain the structure of the peak region, 
which really should be determined by subleading order resummation. To demonstrate the 
danger of trusting a tuned Monte Carlo generator, we run the same event generator at Q — 1 
TeV, and compare again to the theoretical calculations, see Figure [121 Now pythia looks 
like the leading-order event generator that it is, and the hadronization corrections are small, 
but pythia undershoots the more accurate 4th order theoretical prediction. At high energy 
the difference will be more difficult to absorb into non-perturbative effects since hadronization 
corrections are small. One consequence is that these Monte Carlo generators may be under- 
estimating backgrounds at an ILC by 30%, and perhaps by a similar magnitude at the LHC 
as well. 

An alternative to correcting the theoretical distribution with a Monte-Carlo transfer matrix 
is to include explicitly a theoretical model of non-perturbative corrections and then use data 
to determine its parameters. The non-perturbative effects are suppressed by the center-of- 
mass energy and will scale as a power of Anp/Q, with Anp ~ 1 GeV a scale characteristic of 
strong-interaction effects. The effective theory analysis shows that since scales lower than Q 
appear in the perturbative expansion, there will in fact be power corrections suppressed by 
the lowest scale, in this case the soft scale tQ which will go as a power of A NP /(rQ). For 
completely inclusive processes first order power corrections are absent, but one should not 
expect the leading power to be absent for thrust. 

The non-perturbative effects will be most important in the soft region for small r. The 
corrections can be parameterized by a non-perturbative shape function which is convoluted 
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Figure 11: Contours at 95% confidence level for a fit to the opal data of a s and a non- 
perturbative shift parameter A NP . 

with the perturbative soft function [4*8l |4"9] 



Then one can parametrize Snp(&) with a few-parameter family of distributions [50]. For 
example, a common model is £np(&) = S(k — Anp), which leads to an overall shift in the thrust 
distribution. Figure [TT1 shows the result of a simultaneous fit to A NP and a s for the opal data. 
From this rough analysis one can see that the fit to lep data has trouble distinguishing the 
effect of raising the shift parameter from increasing the coupling - both variations increase 
the theoretical prediction in all bins where the fit to data is performed. 

Much of the evidence for a shift in event shape distributions has come from comparisons 
to data of calculations done at NLO or with resummation at NLL (51]. It would be very 
interesting to reconsider these analyses including information from NNLO and with N 3 LL 
resummation. To extract the soft shape function a detailed analysis, including lower energy 
data, should be performed. At lower energies the effect of the power corrections will be more 
pronounced so that the parameters of the shape function can be determined and then used 
in the extraction of a s from higher energy data. The high statistics jade data with energies 
from Q = 22 — 44 GeV might be particularly suitable for such an analysis [52J. 

In our Monte-Carlo studies, we find that quark-mass effects at lep 1 are of order 1%. 
They tend to increase a s , while hadronization effects lower the central value. In fixed-order 
perturbation theory, the quark-mass effects have been evaluated at NLO [531 EH [55J, [56] . Using 
the factorization theorem for the production of massive quark jets [19] and the recent two-loop 
result for the massive jet-function [57J, it is would be possible to perform the resummation 
also for the 6-quark contribution. Since the quark mass corrections are small in the region 
where we extract a s , a fixed-order treatment might be sufficient. Additional issues involved in 
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Figure 12: Comparison between theoretical predictions in effective field theory at first order 
and fourth order, as defined in table HJ and pythia at the parton and hadron level, aleph 
data is included in the first panel. 



matching the perturbative soft and non-perturbative shape functions were discussed recently 
in [58]. 

Since neither Monte-Carlo hadronization corrections nor a simple non-perturbative shift 
model are satisfactory, we conclude that the best option at this point is to fit the parton- 
level distribution. To estimate the hadronization uncertainties, we simply lift the errors from 
previous studies of the aleph and opal data. Numerically this is essentially equivalent to 
using ariadne to calculate the hadronization and quark-mass corrections and the difference 
to pythia as an estimate of the resulting uncertainty, as can be seen in Tables |2] and [31 With 
the increased perturbative precision of our result, it would be important to get better control 
over hadronization effects and to have a more reliable way to assess the associated uncertainty. 
As we discussed above, this can be achieved with a dedicated shape-function analysis involving 
also lower energy data. 



6 Conclusions 

We have resummed the leading logarithmic corrections to the thrust distribution to N 3 LL. 
Our calculation is based on an all-order factorization theorem for the thrust distribution in 
the two-jet region T — > 1. The traditional method for resummation of event shapes is limited 
to NLL. The present paper goes beyond this not only by one but by two orders in logarithmic 
accuracy. 

The factorization theorem, obtained using Soft-Collinear Effective Theory, separates the 
contributions associated with different energy scales in a transparent way. Those associated 
with higher energy scales are absorbed into Wilson coefficients. Solving the renormalization- 
group equations resums large perturbative logarithms of scale ratios. An advantage of the 
effective theory treatment is that the factorization theorem is derived at the operator level. The 
different building blocks in the factorization theorem are given by operator matrix elements 
and appear in a variety of other processes. With the exception of the two-loop constant in the 
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soft function, all the necessary ingredients to the factorization theorem were known to N 3 LL 
accuracy from resummations of other processes. We have determined the missing two-loop 
constant numerically using effective field theory and an NLO fixed-order event generator. 

Comparing to fixed-order results, we found that the logarithmically enhanced pieces, de- 
termined by a few constants in the effective theory, amount to the bulk of the fixed-order 
results, even away from the endpoint T — > 1. Of particular interest is the comparison at 
NNLO. The necessary fixed-order calculation has been completed only recently and so far 
not been independently checked. The close agreement with the logarithmic contributions we 
derive provides a non-trivial check on both calculations. Once matched to the full fixed-order 
result, our result is valid not only to N 3 LL accuracy, but also to NNLO in fixed-order pertur- 
bation theory. Matching improves our result away from the endpoint region, but numerically 
the matching corrections are small, in particular at NNLO. 

Our result is the most precise calculation of an event shape to date, and we have used it 
to perform a precision determination of a s using aleph and opal data. Our final combined 
result is 

a s (m z ) = 0.1172 ± O.OOlO(stat) ± 0.0008(sys) ± 0.0012(had) ± 0.0012(pert) 
= 0.1172 ±0.0022. 

Unfortunately, we had to combine different data sets with the conservative assumption that 
systematic errors are completely correlated. We also had to use the fit regions selected by 
the experiments which are not optimized to our higher order calculation. An improved error 
analysis would involve information from the collaborations about correlations which is not 
publicly available. 

With the resummed calculation, the perturbative uncertainty is finally smaller than the 
other uncertainties at each energy, in contrast to earlier results where the perturbative error 
was dominant. With the reduction in perturbative uncertainty, the hadronization error has 
become a relatively large contribution to the total error. To reduce it, one could parameterize 
the non-perturbative effects with a shape function, and then extract this shape function from 
data at lep and lower energy experiments. In addition, one could account explicitly for quark 
mass effects which should help reduce the systematic errors. 

Even though the perturbative error is greatly reduced by including resummation, the 
technique used to estimate this error may be unduly conservative. We have followed the 
standard procedure and used a collection of scale variations to estimate terms higher order in 
a s . An alternative method, which we have suggested here, is to attempt a more sophisticated 
guess at the effects that a higher order calculation might have. At one higher order the 
resummed distribution is known up to a handful of numbers, such as higher-loop anomalous 
dimensions. So we can extrapolate an approximation to these numbers and use that directly. 
This procedure results in smaller and perhaps more realistic perturbative errors, although we 
have not used the errors derived this way for the final results. 

The effective theory can also be used to study other event shapes. For example, heavy- 
and light-jet masses can be obtained with minimal modifications from the formulae given 
here [151 HH]- These observables involve the same hard and jet functions as (J4j) and the 
necessary soft function can be determined in the same way as we did here for thrust. The 
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factorization theorem for a wider class of event shapes, including jet-broadening and the C- 
parameter was derived recently |59j. Its form is the same as (j3j), except that it involves different 
jet-functions which depend on the variable under consideration. To reach the same accuracy 
we have achieved here, additional perturbative calculations will thus be necessary. 

We could also try to use the same techniques to calculate precision observables in a hadronic 
environment. Many of the necessary ingredients have already been understood from the 
threshold resummation for inclusive processes such as deep-inelastic scattering and Drell- 
Yan production. Despite the complication of hadronic initial states, a precision calculation 
of jet-observables relevant for the LHC seems feasible. Considering the discrepancy we found 
between pythia and the fourth-order effective theory prediction for the thrust distribution 
at 1 TeV (see Figure fT2l) . having a systematically improvable way to perform resummations 
might be vital for the LHC. In addition, given the size of the logarithmic corrections found 
here, it is likely that many fixed-order calculations can be improved using methods of effective 
field theory. 
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A Anomalous dimensions 

The QCD beta-function satisfies 
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where we have written /3 and (3\ in terms of the Casimir invariants Cp 
Tp = §, but have evaluated /3 2 and /3 3 for iV = 3 colors. The RG equation 
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It is also useful sometimes to work with perturbative expansion of ct s (/j,) in terms of a s at a 
fixed renormalization scale, 
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The exact anomalous dimensions are known to a 3 . The anomalous dimensions 
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For the soft function 
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And for the cusp anomalous dimension 
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Analytical expressions for the three-loop terms 7^, 72 and T 2 can be found in [24]. The part 
of the cusp anomalous dimension is not known and we estimate it using a Pade approximation. 
The same approximation works well at and in any case our results are very insensitive to 
the value of r 3 . 



B Hard, jet and soft function 

The hard function can be written as 
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The three-loop constant qf is not yet known but only contributes to the 6(t) part of the thrust 
distribution. The values of the lower order constants are 
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The expression for H(Q, //) is obtained by solving the RG-equation (jSJ) order by order in a s . 
The RG equations for the Laplace transformed jet function and soft function have the same 
form so that their explicit forms are obtained from the above result using simple substitution 
rules. Defining the Laplace transform of the cross section as 
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the cross section factors into the product of the Laplace transforms of the jet- and soft func- 
tions: 
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The Laplace transforms j and St of the jet and soft functions are defined as in ( 1571) . After 
writing these as functions of a logarithm of the argument, the RG equations simplify to 
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Comparing to the RG equation for the hard function Eq. ([5}, and looking at Eq. (154"]) one sees 
that the expression for the jet-function j(L,fi) is obtained from ([55]) . by simple substitutions: 



j(L,n) = h(L, //) with 7 



-7 J , 



and T 



cusp 



-r 



cusp; 



s T (L,n) = h(2L,(i) with j H -> -7 s , and c H 
The constants for the jet and soft functions are 
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The constants cf CF ,cf CA and cf n ^ are extracted numerically as explained in Section [3j We 
found 
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C Singular terms in the thrust distribution 

The fixed-order distributions can be written in terms of delta functions and plus distributions. 

D(r) = 5(r) D s (r) + (g) [D A (r)] + + (^f [D B (r)] + + (^f [D c (r)} + . (65) 
The delta-function terms are known to a 2 s accuracy 
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Table 5: Numerical values for the expansion coefficients of R(t) as defined in (|68|) . 
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The numerical values of cf c , cf c and cf n ^ were given in (138]) . 

To compare with the existing literature and for the readers convenience, we also quote the 
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third order result for the quantity R(t), which is is often written in the form 
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We normalize here to the total hadronic cross section ahad given in (1361) instead of the Born 
cross section uq. Our result provides the normalization of R{t) to second order 
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The numerical values of the above coefficients are listed in Table EJ The NLL coefficients up to 
0(oP s ) were given in [12] and we completely agree with their results. In the same reference the 
remaining a 2 s coefficients were determined using a fit to the numerical fixed order result with 
the result C\ = 34 ± 22 and G21 = 30 ± 10. Our analytical result agrees with the extracted 
value of G21, but our value of G\ is about a factor of two larger. This disagreement is perhaps 
not that surprising, given that [12] had to extract C 2 and G 2 i numerically using a simultaneous 
fit to both quantities at small r, where the result is dominated by the contribution from the 
logarithmic term proportional to G^i- Since we have the analytical result for G21, we are able 
to extract C2 with much higher precision. 
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